Quantum Monte Carlo Study of a Dynamic Hubbard Model 
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The 'dynamic' Hubbard Hamiltonian describes interacting fermions on a lattice whose on-site 
repulsion is modulated by a coupling to a fluctuating bosonic field. We investigate one such model, 
introduced by Hirsch, using the determinant Quantum Monte Carlo method. Our key result is 
that the extended s-wave pairing vertex, repulsive in the usual static Hubbard model, becomes 
attractive as the coupling to the fluctuating Bose field increases. The sign problem prevents us from 
exploring a low enough temperature to see if a superconducting transition occurs. We also observe 
a stabilization of antiferromagnetic correlations and the Mott gap near half-filling, and a near linear 
behavior of the energy as a function of particle density which indicates a tendency toward phase 
separation. 

PACS numbers: 71.10.Fd, 71.30.+h, 02.70.Uu 
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I. INTRODUCTION 

The fermion Hubbard Hamiltonian^, originally pro- 
posed to describe the physics of transition metal monox- 
ides FcO, MnO, and CoO, has been widely used as a 
model of cuprate superconductors, whose undoped par- 
ent compounds, like La2Cu04, are also antiferromag- 
netic and insulating. Indeed, early Quantum Monte 
Carlo (QMC) simulations of the Hubbard Hamiltonian 
suggested that d-wave pairing was the dominant super- 
conducting instability^!^, a symmetry which was subse- 
quently observed in the cuprates^. However, the sign 
problem precluded any definitive statement about a 
phase transition to a d-wave superconducting phase^i^. 
Over the last several years, QMC studies within dynam- 
ical mean field theory and its cluster generalization^ 
are presenting a more compelling case for this tran- 
sition. The existence of charge inhomogeneities in 
Hartree-Fock^ and density matrix renormalization group 
treatments^, along with the experimental observation of 
such patterns^, offer further indications that significant 
aspects of the qualitative physics of the cuprates might 
be contained in the Hubbard Hamiltonian. 

Nevertheless, there are a number of features of high 
temperature superconductors which do not completely fit 
within the framework of the single band Hubbard Hamil- 
tonian. For example, the cuprate gap is set by the charge 
transfer energy separating the copper d and oxygen p 
orbitalsiiii^ as opposed to a Mott gap between copper 
d states split by the on-site repulsion. Considerable ev- 
idence for the possible important role of phonon modes 
in aspects of the physics is availabloi^. 

Hirsch has emphasized the asymmetry in transition 
temperatures, and other properties, between the elec- 
tron and hole doped cuprates as a reason to consider 
more general models, since the particle-hole symmetry 
of the single band Hubbard Hamiltonian requires that 
its behavior be rigorously identical for fillings p = I — x 
and p = I + X. Partially motivated by this asymmetry, 
he introducedi^ii^ii^iiii^ the dynamic Hubbard Hamil- 



tonian, 



H = -tY^ (CjvCi<T + cLcj^) - fJ-^ini] + nil) 

(iJ> i 

+ X] [^o<yi + a^oCTi +{U - 2gujoai ) ni^nn .(1) 

i 

Here the first term, involving the fermion creation (de- 
struction) operators Cj^(cjo-) at site j with spin a, is the 
tight binding kinetic energy describing the hopping of 
electrons between near neighbor sites. We consider here 
a two-dimensional square lattice and chose t = I to set 
our scale of energy. The on-site interaction energy dif- 
fers from the usual static Hubbard Hamiltonian in that 
its value U is modulated by a dynamic field cr? which 
can take the values a? — ±1. As a consequence the on- 
site repulsion U has bimodal values t/min = U — Sgwo 
and C/niax = U -\- 2guj{j. This dynamic field itself has 
non-trivial quantum fluctuations controlled by the rela- 
tive values of the longitudinal and transverse frequencies 
gujQ and lo^). Here cr? and af are Pauli matricesi^. The 
variation in U , Hirsch argued, has its physical origin in 
the relaxation which occurs with multiple occupation of 
an atomic orbital. 

Hirsch and collaborators have studied the physics of 
Eq.[T]with a variety of methods, including a Lang-Firsov 
transformation (LFT)^^, exact diagonalization (ED) of 
small clustersi^iii, and world-line Quantum Monte Carlo 
(WLQMC)^ in one dimensioi^. Within the LFT it is 
seen that the hopping of electrons is renormalized by the 
overlap of the states of the dynamic variable on neighbor- 
ing lattice sites. Superconductivity then arises because 
isolated holes are essentially localized by a small overlap, 
whereas holes that are on the same or neighboring sites 
can move around the lattice. Furthermore this effect is 
operative for holes in a nearly flUed system, but not elec- 
trons in a nearly empty lattice. Thus pairing is linked to 
the presence of holes, and the physics is manifestly not 
particle-hole symmetric. ED provided quantitative val- 
ues for the overlaps and confirmed the picture based on 
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FIG. 1: (Color online) Binding energy UbB on a four site 
cluster as a function of the coupling g to the dynamically 
fluctuating field . The binding energy f/cff can go negative 
at large g, suggesting the possibility of pairing. 

the LFT on small clusters. 

ED also allows for the evaluation of the 'binding en- 
ergy', C/eft = 2Eo{N + 1) - Eo{N + 2) - £;o(^). Here 
Eo{N) is the ground state energy of a cluster with N 
electrons. A negative UeS indicates that it is energeti- 
cally favorable to put two particles together on a single 
cluster rather than separate them on two different clus- 
ters. On a sufficiently large lattice, two particles would 
tend to be close spatially rather than widely separated. 
In Fig. [1] we show an evaluation of [/off on a 2x2 lattice. 
These numbers were obtained independently from, but 
arc identical to, those of Ref. [H. As the coupling g to 
the dynamic field increases, UeS is driven negative, indi- 
cating the possibility of binding of particles and hence su- 
perconductivity. WLQMC simulations in one dimension 
confirmed this real space pairing by explicitly showing 
the preference of the world lines of holes to propagate 
next to each other and a large gain in kinetic energy 
when the hole-hole separation becomes small. Signifi- 
cantly, these simulations also showed that the kinetic en- 
ergy disfavors proximity of holes in the Holstcin model, 
which also features the tendency of holes to clump to- 
gether by distorting a local phonon degree of freedom. 
Thus pairing in the dynamic Hubbard model is distin- 
guished from that of more traditional electron-phonon 
models by being driven by the kinetic energy as opposed 
to a potential energy. 

In this paper we examine the properties of the dynamic 
Hubbard Hamiltonian with determinant Quantum Monte 
Carlo (DQMC)2i. This approach allows us to work in 
two dimensions, as opposed to previous {d — 1) WLQMC 
studies, and also to examine lattices of an order of magni- 
tude greater number of sites than ED. On the other hand, 
the ability of DQMC to reach low temperatures is limited 
by the sign problen*^. We find that the extended s-wave 
pairing vertex, which is repulsive in the static Hubbard 
model, is attractive in the dynamic model, that is, ex- 
tended s-wave superconducting correlations are enhanced 
by the dynamic fluctuations. However, the pairing sus- 



ceptibilities are still only rather weakly increasing down 
to the lowest temperatures accessible to us (temperature 
T greater than 1/40 the electronic bandwidth). 

We also find, near half-filling, that the antiferromag- 
netic correlations can be enhanced relative to the static 
Hubbard Hamiltonian, particularly for densities above 
p = 1. The Mott gap can also be stabilized. Interest- 
ingly, the total energy appears to be close to linear in the 
particle density, as opposed to a clear concave up curva- 
ture in the static Hubbard model (with either repulsive 
or attractive interactions). 

The organization of this paper is as follows: In the next 
section we present our computational method, DQMC, as 
it applies to the dynamic Hubbard model. We describe 
several minor adjustments to the DQMC algorithm for 
the static Hubbard model that are needed in order to 
study the dynamic model. Our observables arc also de- 
fined. In Section HI, we present the results from our 
Monte Carlo simulations. The topics of antiferromag- 
netism and the Mott transition, pair susceptibilities and 
superconductivity, and the energy characteristics of the 
dynamic Hubbard model are discussed. The paper closes 
with conclusions in Section IV. 



II. COMPUTATIONAL METHODS 

Although he did not undertake such studies, Hirsch 
pointed ouli^^ that the dynamic Hubbard model could 
be simulated with a relatively minor modification of the 
DQMC method^!. In DQMC, an auxihary 'Hubbard- 
Stratonovich' (HS) field is introduced to decouple the 
on-site Hubbard repulsion. The trace over the resulting 
quadratic form of fermion operators is performed ana- 
lytically, leaving an expression for the partition function 
which is a sum over the HS variables whose weight is 
given by the product of two determinants, one for spin 
up and one for spin down, that are produced by evaluat- 
ing the trace. 

In DQMC for the usual Hubbard Hamiltonian, the HS 
field couples to the difference between the up and down 
spin electron densities, with a coupling constant which is 
independent of spatial site and imaginary time. In a sim- 
ulation of the dynamic Hubbard model, the coupling of 
the HS field depends on the dynamic field cr?{T). The 
imaginary time dependence arises from the transverse 
term af in the Hamiltonian. When the path integral for 
the partition function is constructed cr? induces flips be- 
tween the two values a? = ±1, so this quantity becomes 
dependent on r. As a consequence, minor modifications 
are required to the standard expressions2i for the ratio 
of determinants before and after the Monte Carlo move 
of a HS variable, and for the re-evaluation of the Green's 
function. 

The dynamic field variables must also be updated, and 
again, only minor modifications of the formulae for the 
determinant ratio and Green's function update are re- 
quired. A final difference is that there is a contribution 
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FIG. 2: (Color online) Double occupancy (left) and expecta- 
tion value of dynamic field (right) as functions of the cou- 
pling g. The solid line is the result of exact diagonaliza- 
tion and the symbols of the determinant QMC simulations. 
The cluster size is 2x2 (the same as for the binding en- 
ergy calculation of Fig. [T] and Ref. [ish . Parameters are 
t = 1, [/ = 4, /? = 1.3, = 2 and cjo = 1.5. 

to the weight coming from the and terms in the 
Hamiltonian. The former try to ahgn the dynamic vari- 
ables in the imaginary time direction, while the latter 
favor positive values of the dynamic field. Such pieces 
of the action, which enter the weight of the configuration 
along with the fermion determinants, are similar to those 
arising in simulations of the Holstein Hamiltonian^. 

We verified our DQMC code by comparing to exact di- 
agonalization results on a 2x2 spatial lattice (Figs. [TT^ . 
and also by checking analytically soluble limits such as 
t ~ 0. The results of our DQMC/diagonalization cal- 
culations on 2x2 lattices are completely consistent with 
those of Hirsch. For example, we have quantitatively re- 
produced the binding energy plot. Fig. l(top) of Ref. [Tsl 
and our Fig. [1] As a further check, we compared DQMC 
results for the double occupancy, (n-^ni), and the expec- 
tation value of the dynamic field, (cr^), to results from 
ED. See Fig. [2 

We did not observe any major difference in the charac- 
teristics of the DQMC algorithm in simulating the dy- 
namic Hubbard model: Autocorrelation times remain 
short, as is typically the case with DQMC, and there 
was no major change in the numerical stability^ i^'^i^^i^^ . 
The key issue in DQMC is the 'sign problem' which we 
will discuss in the following sections. 

DQMC allows us to measure any observable which can 
be expressed as an expectation value of products of cre- 
ation and destruction operators. Our measurements in- 
clude the energy (H) (not including the chemical poten- 
tial term), particle density p = (n), and Green's function 
Gij(r) = (ci(r)cj(O)), as well as the average of the dy- 
namic field (erf). The dependence of the density on the 
chemical potential fj, and the Green's function, when an- 
alytically continued to the spectral function, allows us to 
examine, among other things, the Mott metal-insulator 



transition. 

In addition to these single particle properties we also 
examine magnetic correlations, and specifically, the mag- 
netic structure factor, 

S{k) - J2 ( ("J+'T - '^j+u)("jT - ) ■ (2) 
1 

Our focus will be on the antiferromagnetic response, 
5(k= (tt,^)). 

We look at superconductivity by computing the corre- 
lated pair field susceptibility. Pa, in different symmetry 
channels. 

Pa = / dr(A„(r)Al(0)) 
Jo 

k 

/.(k) = 1 

/s*(k) = cos{kx) + cofi{ky) 
fd{k) = cos{kx) - cos(fcy) . (3) 

These quantities can also be expressed in real space, 

A^ = J^I14^{4+H+4+yi+4~H+4-y^ 

A^ = ^ E 4, {4^H 4^yi + - -(4) 

The correlated susceptibility P^ takes the expectation 
value of the product of the four fermion operators enter- 
ing Eq. [3l We also define the uncorrelated pair field sus- 
ceptibility Pa which instead computes the expectation 
values of pairs of operators prior to taking the product. 
Thus, for example, in the s-wave channel, 

= jpT. f dT{c,,{T)c,,{T)c\,{Q)cl{Q)) 

- 4E f dT{c,,{r)cl{Q)) (ciT(r)cyO)) (5) 

Pa includes both the renormalization of the propagation 
of the individual fermions as well as the interaction ver- 
tex between them, whereas Pa includes only the former 
effect. Indeed by evaluating both P and P we are able 
to extract^ the interaction vertex F, 



If Fq, < 0, the associated pairing interaction is attractive. 
Fq — 1 signals a superconducting instability. 
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III. RESULTS 

A. Mott Transition and Antiferromagnetism 

It is useful to begin our study of the dynamic Hubbard 
model by understanding the behavior of the dynamic field 
at different fillings (Fig. [3|) . For fillings below one particle 
per site, p < 1, the dynamic field af ~ — 1 because of the 
coupling to the external ficldgwo hence the interaction 
U + 2guJo ~ i7max and the double occupancy is reduced. 
However, once double occupancy is unavoidable (p > 1), 
the interaction term strongly favors erf = +1. Fig. [3] 
shows that this evolution from negative to positive values 
is nearly linear once p > 1. Meanwhile, the expectation 
value of erf measures the fluctuations of erf in imaginary 
time. It is not surprising, then, that this quantity ex- 
hibits a maximum at roughly the midpoint between the 
evolution from erf = —I to erf = +1, at p w 1.5. In 
the results of Fig. [31 and throughout this paper unless 
otherwise stated, the simulations were performed on 6x6 
lattices. 

Next, we compare the Mott gap and magnetic corre- 
lations in the static and dynamic Hubbard models. In 
Fig. [?] we plot the density p as a function of chemical 
potential p. A plateau at p = 1 indicates the formation 
of a Mott insulator. The cost to add a particle sud- 
denly jumps by U because additional particles are forced 
to sit on sites which are already occupied. At the in- 
verse temperature chosen, /3 = 5, for the static Hubbard 
model, the plateau is only beginning to develop. How- 
ever for the dynamic model the plateau is much more 
robust. This is expected since near half- filling, as we 
have seen, the on-site repulsion mostly takes on its max- 
imum value i/inax = 7.8, for the parameters in Fig. |4l 
We have chosen dynamic Hubbard parameters g and ojq 
which gets the system as close as possible to the most at- 
tractive (negative) binding energy [/off while still having 

C/min > 0. 

Figure [S] gives further insight into the behavior of the 
density near full filling. In the static model, the cost to 
add particles to the system is set by the on-site U (in 
the case that U exceeds the bandwidth W — ^t). How- 
ever, in the dynamic model, as full filling is approached, 
the double occupancy cost is reduced to C/min- For the 
parameters chosen in Fig. O C/min = 0.2 is close to zero. 
Thus we expect the filling of the lattice to be complete 
when the chemical potential reaches the top of the band. 
At, in good agreement with the plot. 

The static Hubbard model exhibits antifcrromagnctic 
correlations at half-filling on a bipartite lattice, since only 
electrons with anti-aligncd spins can hop between neigh- 
boring sites. This leads to a lowering of the energy by the 
exchange energy J = 4t^/C/ relative to sites with parallel 
spin, for which hopping is forbidden. Indeed, a finite size 
scaling analysis of the structure factor has shown there 
is long range order in the ground state^l. Figure [B] com- 
pares the value of the antifcromagnetic structure factor 
S{tt, tt) for the static and dynamic models. At half-filling, 




FIG. 3: (Color online) (cr'^) and (cr"^) as a function of p for 
ujo = 0.5, g = 3.8. From p = to half-filling, the system 
minimizes its energy by maximizing the interaction term U — 
2gLuoa^ to avoid double occupation, that is, {a^) ~ —1. In 
this Figure, and elsewhere in this paper, the lattice size is 6x6 
unless otherwise stated. 
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FIG. 4: (Color online) Comparison of the evolution of the den- 
sity p with chemical potential p for the static and dynamic 
Hubbard models. The dynamic model has a significantly bet- 
ter developed Mott insulating gap, as well as a pronounced 
particle-hole asymmetry. Here the inverse temperature f3 = 5. 



S(tt, tt) for the dynamic model attains a maximal value 
50% larger than that of the static model. There is a 
marked asymmetry in the magnetic response at values 
greater and lower than p = 1 in the dynamic model, with 
S{'!T, tt) remaining high to values of p ten percent larger 
than half-filling. We also show results for the negative 
U Hubbard model, which has no tendency for magnetic 
order at any filling. (Instead, the attractive Hubbard 
model exhibits long range charge and superconducting 
correlations at p = 1). 

The spectral function A{uj), which we obtain with an 
analytic continuation of G{t) using the maximum en- 
tropy method^, shows supporting evidence for the en- 
hancement of the Mott gap at half-filling. Fig. [T^top). 
Above half-filling A{uj) exhibits a sharp resonance at 
ti> = 0, Fig. [TKbottom). The comparison of A{uj) for 
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FIG. 5: (Color online) The density p as a function of chemical 
potential p at [/ = 4 and /? = 5. As the coupling g increases, 
the cost to add particles to an already occupied site decreases. 
As a consequence, p rises more steeply with ^. 
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FIG. 6: (Color online) The antiferromagnetic structure factor 
S{tv, tt) at inverse temperature /9 = 5 as a function of density 
p. For both the static and dynamic repulsive Hubbard Hamil- 
tonians there is significant antiferromagnetic order near half 
filling, with the magnetic correlations in the dynamic model 
somewhat more robust. There is no magnetic signal for the at- 
tractive model which, instead, is known to show strong charge 
density wave and s-wave superconducting correlations. 



p — 0.5 and p ~ 1.5 further emphasizes the lack of 
particle- hole symmetry, Fig. [Tfbottom) . 
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FIG. 7: (Color online) Top: Comparison of the spectral func- 
tion A{lo) for the static and dynamic Hubbard models at 
/3 = 5 and half-filling. We can see clearly that the Mott 
gap is more robust in the dynamic case. Middle and bottom: 
The behavior of A{lo) away from half-filling22.. In all cases 
the spectral function is finite at the Fermi energy, indicating 
metallic behavior. However, for the dynamic model a,t p — 1.5 
there is a sharp resonance aX uj = whereas in the other cases 
the spectral function is suppressed there. 



B. Pairing Susceptibilities 

We turn now to a discussion of superconductivity in 
the dynamic model. In the static Hubbard model, it has 
been shown that the s-wave pairing vertex is repulsive 
(positive). The d-wave vertex is negative, but only rel- 
atively weakly so at the temperatures accessible to the 
simulationa^i^. Near half-filling, the extended s-wave ver- 
tex is also attractive, but markedly less so than c?-wave, 
suggesting that d-wave symmetry is the most likely in- 
stability. However, the same sign problem which pre- 



cludes a definitive statement about superconductivity in 
the static model also limits what we can conclude here for 
the dynamic model. Nevertheless, there is an interesting 
qualitative difference between the two models which can 
be clearly discerned. 

Specifically, the extended s-wave vertex is attractive 
in the dynamic model in the regime of g where UcS 
is negative, while it is repulsive in the static model at 
these high fillings. In Fig. [5] we compare the temper- 
ature evolution of the correlated and uncorrelated sus- 
ceptibilities, Ps' and Ps*, at p = 1.89 and U = 3 and 
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FIG. 8: (Color online) The extended s-wave pair susceptibili- 
ties Ps and Ps* as a function of temperature for U = 3,ijJo ~ 
0.5, and g — 2.9. Here, unlike the static model, Ps* exceeds 
Pa* when the temperature is lowered. However, we cannot 
say if Ps* might diverge at low temperature because of the 
sign problem. 
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FIG. 9: (Color online) Ps* as a function of density p for (3 = 5 
near full-filling. The s* channel becomes attractive when g 
increases. Fs* ^ — 1 would signal a superconducting instabil- 
ity. 



see that the Pg* > Ps* ■ The average sign takes the val- 
ues 0.94, 0.92, 0.83, 0.73, and 0.63 at /3 = 5,6,7,8,9 
respectively. The resulting attractive (negative) vertex 
is given in Fig. O For g = 0, the static model, the ver- 
tex is repulsive. But it systematically decreases and goes 
negative as the coupling to the dynamic field is strength- 
ened. In this plot the inverse temperature is fixed at 
/? = 5 and the density is allowed to vary. The average 
sign takes the values 0.36, 0.54, 0.73, 0.89, and 0.96 at 
p= 1.50, 1.60, 1.70, 1.80, L90. 

Figure [TO] (top) shows that, in contrast to the behavior 
of Ts* , the s-wave vertex is strongly repulsive, although 
g does weaken the repulsion somewhat as it increases. 
Meanwhile, we see in FiglTO] (bottom) that near full filling 
the d-wave vertex is more weakly attractive than the s*- 
wave vertex. This suggests that if the dynamic Hubbard 
model does have a superconducting instability at small 
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FIG. 10: (Color online) Top: Fs as a function of p for different 
values of g at 13 — 5. Unlike Fs* , the s-wave channel remains 
repulsive. Bottom: F^j ^ ^ as a. function of p. The d^2_y2- 
wave channel is attractive, but the effect is less pronounced 
than for s*, especially near full filling. 



hole-doping that it would be of s* symmetry, unlike the 
d-wave symmetry which is most attractive for the static 
model^. 

It is informative to compare the onset of attraction 
in the pairing vertex with the development of negative 
binding energy. Fig. Illl shows Fc. T^* and F^ ^ ^ versus 
5 for [/ = 4 and ujq = 0.5. The filling p = 1.8." On 2x2 
lattices, for which the ED calculation of C/eff is feasible, 
Fs« becomes negative at somewhat larger values of g than 
where C/cff becomes negative. The figure also shows that 
F is relatively insensitive to lattice size: the 2x2 and 
6x6 lattices give results which arc quantitatively rather 
similar for most values of g. Note also that F^* is strongly 
repulsive in the static model g = 0. 

A significantly larger enhancement of superconductiv- 
ity was reported^ in a Hubbard Hamiltonian in which 
the hopping of one spin species is modulated by the den- 
sity of the other. This model was argued to be connected 
to the dynamic Hubbard Hamiltonian in the limit of large 
Wo- We conclude this section by exploring the wq depen- 
dence of the pairing vertex, to see if larger wq might show 
a greater tendency for superconductivity. In Fig. [TO] we 
show the vertices as a function of wq- We have fixed the 
product gwo = 1-9 and C/ = 4 so we can stay near the 
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FIG. 11: (Color online) Pairing vertices and binding energy as 
functions of the dynamical coupling g. Ts* becomes negative 
at t; > 3.5 while C/eff < when g > 2.3. As long as g < 4, 
i.e. for the entire range of the horizontal axis, both C/min and 
Umax are repulsive. T^' is strongly repulsive in the static 
model 3 = 0. 



0.6 



0.4 



0.2 



-0.2 





1 1 1 1 1 1 1 






p=1. 


8 U=4.0, gco„=1.9 1 








T *" 




N=4x4 








N=4x4 








N=4x4 


1 , , , , 





0.5 



1.5 



FIG. 12: (Color online) Pairing vertices as functions of fre- 
quency Wo at fixed gwo = 1.9, p — 1.8, (3 = 5, and U = 4. 
The attractive d-wave vertex shows only a weak dependence 
on wq. 
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FIG. 13: (Color online) Total energy as a function of p at 
/3 = 5 for the static attractive and repulsive Hubbard models, 
and for the dynamic model. The static models show clear 
positive curvature, indicative of thermodynamic stability. In 
the dynamic model E{p) is nearly linear. 
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FIG. 14: (Color online) Total energy as a function of p at 
(3 = 5. For g = 0, the static model, the curvature is positive. 
As g gets larger, the energy becomes linear in the density. 



C. Energy 

The total energy (Fig fT5|) also shows a markedly differ- 
ent dependence on the density, p, in the dynamic Hub- 
bard Hamiltonian. Whereas the static positive and neg- 
ative U Hubbard Hamiltonians have cPE/dp^ > 0, the 
positive curvature of the dynamic model that is evident 
below half-filling becomes very small for p > 1 as g 
increases and eventually the curvature nearly vanishes. 
Figure fT4l shows this linear behavior developing with g. 

The temperatures at which we performed our simula- 
tions are low enough that the total internal energy, E, is 
nearly equal to the free energy, F. As it is well known, 
negative curvature in the free energy as a function of the 
density, in the canonical ensemble, leads to negative com- 



pressibility and is thus a signal for phase separation and a 
first order phase transition^. Thermodynamic stability 
requires positive curvature for the free energy versus den- 
sity. While our simulations are performed in the grand 
canonical ensemble, where such negative curvatures are 
not observed, we do see (Fig. I14p a progression from pos- 
itive to zero curvature as 17 -^■ 4. At the same time, and 
recalling that /i = d{F/V)/dp, we sec in Fig. [5] that as 
(7 — > 4, the p versus p curves get steeper signalling higher 
compressibility k = dp/ dp. Noting that C/min vanishes 
for (7 = 4 and becomes negative when g > 4, we interpret 
these observations as a possible phase separation setting 
in at g = 4 whereby the system develops hole-rich and 
hole-deficient regions. 
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IV. CONCLUSIONS 

111 this paper we have performed determinant Quantum 
Monte Carlo simulations of a two dimensional Hubbard 
Hamiltonian in which the on-site repulsion is coupled to 
a fluctuating bosonic field. Our studies complement ear- 
lier work using the Lang-Firsov transformation and exact 
diagonalization and QMC in one dimension. We note a 
number of interesting features of the model. First, the 
Mott gap at half-filling is stabilized. Second, antifcr- 
romagnetic correlations are enhanced above half-filling. 
The extended s-wave pairing vertex, which is repulsive in 
the ordinary static Hubbard Hamiltonian, is made attrac- 
tive in the dynamic model. The value of g for which this 
attraction manifests is roughly consistent with the value 



at which the binding energy UaS goes negative on 2x2 
clusters. The sign problem prevents simulations at low 
temperatures to see if an actual pairing instability occurs. 
We have also observed that as g — > 4, i.e. as C/min ~^ 0, 
E{p) becomes linear in p signalling possible phase sep- 
aration into regions of holc-dcflcicnt and hole-rich re- 
gions when C/niin becomes negative for g > A. Finally, 
we note that we have also found, within the Hartree- 
Fock framework, that charge inhomogeneities (stripes) 
are supported by this dynamic Hubbard model^. 
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